#实验参考文心一言
# 安装和载入Seurat包
install.packages("Seurat")
install.packages('dplyr','patchwork')
install.packages("devtools")
devtools::install_github("USERNAME/FilterCells")

# 读取数据
library(Seurat)
library(dplyr)


pbmc.data <- Read10X(data.dir = "D:/Desktop/R语言/第七次/filtered_gene_bc_matrices/hg19")
#进行数据读取

pbmc <- pbmc <- CreateSeuratObject(counts = pbmc.data, 
                                   project = "pbmc3k",)
pbmc[['percent.mt']] <- PercentageFeatureSet(pbmc, pattern = "^MT-|^mt-")

# 数据QC
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)


VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), 
        ncol = 3)

# 归一化
pbmc <- NormalizeData(pbmc)

# 找到高变基因
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)

# 聚类分析
pbmc <- ScaleData(pbmc)
pbmc <- RunPCA(pbmc, npcs = 30)
pbmc <- RunUMAP(pbmc, dims = 1:30)
pbmc <- FindNeighbors(pbmc, dims = 1:30)
pbmc <- FindClusters(pbmc, resolution = 0.5)



# 绘制UMAP图####
DimPlot(pbmc, reduction = "umap", group.by = "ident")



# 绘制每个cluster的Top 2基因气泡图####
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, 
                               min.pct = 0.25, logfc.threshold = 0.25)


top_genes <- pbmc.markers %>%  
  group_by(cluster) %>%  
  slice_max(n = 2, order_by = avg_log2FC)  

# 确保您的Seurat对象名为pbmc，并且您想要绘制RNA assay的数据  
# 提取每个cluster的top基因名称  
top_genes_list <- top_genes$gene  
output_dir <- "D:/Desktop/R语言/第七次/results1" # 定义文件夹名  
output_filename_prefix <- "Cluster_DotPlot" # 定义文件名前缀  

# 确保文件夹存在，如果不存在则创建它  
if (!dir.exists(output_dir)) {  
  dir.create(output_dir)  
}  

# 定义完整的文件路径（包括文件名前缀和扩展名）  
output <- file.path(output_dir, paste0(output_filename_prefix, ".pdf")) 
# 绘制DotPlot  
library(ggplot2)
pdf(file = paste0(output, "_Cluster_DotPlot.pdf"), width = 8, height = 6)  
p10 <- DotPlot(pbmc, group.by = "seurat_clusters", features = top_genes_list) +  
  theme(axis.text.x = element_text(angle = 45, hjust = 1))  
print(p10)  
dev.off()

# 挑选几个marker gene并绘制feature plot和VlnPlot#######
# 假设我们挑选了cluster 0的top 1 marker gene "gene1" 和 cluster 1的top 1 marker gene "gene2"
feature_genes <- c("FGFBP2", "GZMB")

# 绘制feature plot
library(ggplot2)
FeaturePlot(pbmc, features = feature_genes, reduction = "umap")

# 绘制VlnPlot####
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), 
        ncol = 3)
VariableFeaturePlot(pbmc)

